#include "Random.hpp"
#include <iostream>
#include <armadillo>
using namespace std;

namespace Resample{

using namespace arma;

inline void cumsum(double* vect,int N)
{
	int i=1;
	for(i=1;i<N;i++)
	{
		vect[i]=vect[i-1]+vect[i];
	}
}

inline void normalize(double* vect,int N)
{
	double sum=0;
	int i=0;
	for(int i=0;i<N;i++)
	{
		sum=sum+vect[i];
	}
	for(int i=0;i<N;i++)
	{
		vect[i]=vect[i]/sum;
	}
}

class SysResample
{
	public:
		void operator()(int *N,double* prob, mat &vect)
		{
			int n=N[0];
			int m=vect.n_rows;
			mat res(m,1);
			int i=0;
			double temp=0;
			int k=0;
			double c=0;
			double u=0;
			double u2=0;
			//GetRNGstate();
			boost::random::uniform_01<> U;
			RandomG::Random<boost::random::uniform_01<> > unif_rand(U);
			u=unif_rand()*(double)1/n;
			normalize(prob,n);
			cumsum(prob,n);
			u2=u;
			while(i<n)
			{
				if(u2<prob[k])
				{
						
					res(i,0)=vect(k,0);
					u2=(double)(i+1)/n+u;
					i++;
				}else{
				k++;
				}
			}
		//	return res;
			vect=res;
		}

};


class Multinomial{
	public:
	// n multinomial drawing from the set {1,...,*n} 
	// prob: n-dimensional vector of probabilities
	// Remark: O(n) algorithm
	void operator()(int *n, double *prob, mat &res)
	{
		int i,k;
		double sim[*n+1];
		mat foo(*n,1);
		
				
		boost::random::exponential_distribution<> U;
		RandomG::Random<boost::random::exponential_distribution<> > exp_rand(U);
		

		for(k=0;k<=*n;k++)
		{
		 sim[k]=exp_rand();
		}
		cumsum(sim,*n+1);


		for(k=0;k<*n;k++)
		{
		    sim[k]=sim[k]/sim[*n];
		}
		normalize(prob,*n);
	    	cumsum(prob,*n);
		i=0;
		k=0;
		while(k<*n && i<*n)
		{

		       while(sim[i]<=prob[k] && i<*n)
		      {
			  foo(i,0)=res[k];
			  i++;
		      }
	    	k++;

		}
		res=foo;

	}
};

	
class Residual{
	public:
		void operator()(int *N, double *prob, mat &vect)
		{
			int R=0;
			int n=N[0];
			double *wbar=new double[n];
			normalize(prob,n);
			for(int i=0;i<n;i++)
			{
				R+=floor(n*prob[i]);
			}
			for(int i=0;i<n;i++)
			{
				wbar[i]=(n*prob[i]-floor(n*prob[i]))/(n-R);
			}
			Multinomial M;
			int nr=n-R;
			mat vect2=vect;
			for(int i=0;i<n;i++)
			{
				wbar[i]=(n*prob[i]-floor(n*prob[i]))/(n-R);
				//cout << "wwwwwwwww "<< wbar[i] << "   ";
			}
			M(&nr,wbar,vect2);
			for(int i=0;i<nr;i++) 
			{
				vect(i,0)=vect2(i,0);
			}
			int add;
			int kk=nr;
			for(int i=0;i<n;i++)
			{
				add=floor(prob[i]*n);		
				for(int j=0;j<add;j++)
				{
					vect(kk,0)=i;
					kk++;
				}
						
			}
	


			delete[] wbar;
		}

};


class condMultinomial{
	public:
	mat operator()(int *N, double *prob, mat z)
	{
		boost::random::uniform_01<> U;
		RandomG::Random<boost::random::uniform_01<> > unif_rand(U);
		int n=N[0];
		int p=z.n_rows;
		normalize(prob,p);
		cumsum(prob,p);
		mat res(n,1);
	 	//cout << z;	
		for(int i=0;i<n;i++)
		{
			double u=unif_rand();
			int k=0;
			while(u>prob[k])
			{
				k++;
			}
			res(i,0)=z(k,0);	
			
		}
		return res;
	}
};
	
}//end namespace

